{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Sampling and Hypothesis Testing\n",
    "\n",
    "## Carbon-14 assignment: suggested solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "np.set_printoptions(precision=2)\n",
    "\n",
    "# To make the working a little clearer, I have not rounded the results presented below.\n",
    "\n",
    "# Of course, numerical results should normally be presented at an appropriate precision.\n",
    "# In this exercise, none of the values given has more than 2 significant figures, so that \n",
    "# would be a sensible precision for the final results for each question.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "You have been given 500g of wood charcoal and a [Geiger-Müller (G-M) counter](https://en.wikipedia.org/wiki/Geiger_counter). \n",
    "\n",
    "You are using this equipment to investigate the decay of the naturally-occurring radioactive isotope, Carbon-14.\n",
    "\n",
    "You point the G-M tube towards the charcoal.\n",
    "\n",
    "Over six consecutive ten-minute intervals, the detector gives the following counts:\n",
    "10, 21, 19, 15, 26, 17.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1\n",
    "\n",
    "Based on these observations, calculate the\n",
    "\n",
    "(a) sample mean,\n",
    "\n",
    "(b) sample standard deviation, and\n",
    "\n",
    "(c) standard error of the mean\n",
    "\n",
    "for the number of counts observed in a 10-minute interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xbar:  18.0\n",
      "std:  5.440588203494177\n",
      "stderr:  2.221110833194358\n"
     ]
    }
   ],
   "source": [
    "data = np.array([10,21,19,15,26,17])\n",
    "\n",
    "xbar = data.mean()\n",
    "print('xbar: ', xbar)\n",
    "\n",
    "# using the unbiased estimator\n",
    "std = data.std(ddof=1)\n",
    "print('std: ', std)\n",
    "\n",
    "n = len(data)\n",
    "stderr = std / np.sqrt(n)\n",
    "print('stderr: ', stderr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2\n",
    "\n",
    "(a) Plot the PMF for the number of counts observed in a ten-minute interval. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x1a21ebc450>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASo0lEQVR4nO3df6xkd1nH8feHu+2CQIq0i4H+cFtawVsh2F4LRERCFVuTsiDFbBHtH00KShMNEmlNDNCEP0qUmmgjqbbQoNg2VXQ3oJVY/BGCtXehUJZ2dfcKdCmhW1uKaKRs+/jHOSvT25l7z+7O3pk5834lNzPnnOfeefbszGfOfM+cc1JVSJL662mTbkCSdGwZ9JLUcwa9JPWcQS9JPWfQS1LPbZp0A6uddNJJtXXr1km3IUkzZdeuXQ9V1ZZhy6Yu6Ldu3cry8vKk25CkmZLkq6OWOXQjST1n0EtSzxn0ktRzBr0k9ZxBL0k9Z9BLUs8Z9FIXKytw9tmwaVNzu7JyZDXSBBj0UhcXXQT33QePP97cXnTRkdVIE2DQS13s2QNPPNHcf+KJZvpIaqQJMOilLl70orWnu9ZIE2DQS13s3Ln2dNcaaQIMeqmLM85Ye7prjTQBBr0k9ZxBL0k9Z9BLUs8Z9JLUcwa9JPWcQS9t5KkLPE2CJsCglzby1AWeJkETYNBLG3nqAk+ToAkw6KWNPHWBp0nQBBj00kaeusDTJGgCDHppI09d4GkSNAEGvST1nEEvST1n0EtSzxn0ktRzBr0k9ZxBL0k9Z9BLUs8Z9JLUcwa9JPWcQS9JPWfQS1LPdQr6JBck2ZNkb5IrhyzfnOSWdvmdSba2849LclOSe5Lcm+Sq8bYvSVrPukGfZAG4DrgQWAQuSbK4quwy4JGqOhO4Frimnf9mYHNVvQQ4F3jboTcBSSN4FSqNWZct+vOAvVW1UlWPATcD21bVbANuau/fBpyfJEABz0yyCXgG8Bjw7bF0LvWVV6HSmHUJ+pOB+wem97fzhtZU1UHgUeBEmtD/b+AbwNeA362qh1c/QJLLkywnWT5w4MBh/yOkXvEqVBqzLkGfIfOqY815wOPAC4DTgd9M8pQTcFfV9VW1VFVLW7Zs6dCS1NEsDoN4FSqNWZeg3w+cOjB9CvDAqJp2mOYE4GHgLcDfVtX3qupB4DPA0tE2LXU2i8MgXoVKY9Yl6O8CzkpyepLjge3AjlU1O4BL2/sXA3dUVdEM17w2jWcCrwDuG0/rUgezOAziVag0ZusGfTvmfgVwO3AvcGtV7U5ydZLXt2U3ACcm2Qu8Ezj0FczrgGcBX6J5w/hwVX1xzP8GaTSHQSTSbHhPj6WlpVpeXp50G+qLlRV44Qu/P71v3/At5AzsZhr1mpi2GmlAkl1VNXRo3CNj1W8Og0gGvST1nUEvST1n0EtSzxn0ktRzBr0k9ZxBL0k9Z9BLUs8Z9JLUcwa9JPWcQS9JPWfQS1LPGfSS1HMGvST1nEEvST1n0EtSzxn00iyaxYuea2IMemkWzeJFzzUxBr00i2bxoueaGINemkVe9FyHwaCXZtHOnWtPSwMMemkWedFzHQaDXpJ6zqDX7PIrhlInBr1ml18xlDox6DW7/Iqh1IlBr9nlVwylTgx6zS6/Yih1YtBrdvkVQ6kTg16Ses6gl6SeM+glqecMeknquU5Bn+SCJHuS7E1y5ZDlm5Pc0i6/M8nWgWUvTfLZJLuT3JPk6eNrX5K0nnWDPskCcB1wIbAIXJJkcVXZZcAjVXUmcC1wTfu7m4A/Bd5eVWcDrwG+N7buJUnr6rJFfx6wt6pWquox4GZg26qabcBN7f3bgPOTBHgd8MWq+gJAVf1nVT0+ntYlSV10CfqTgfsHpve384bWVNVB4FHgROBHgEpye5LPJfmtYQ+Q5PIky0mWDxw4cLj/BknSGroEfYbMq441m4BXAb/U3r4xyflPKay6vqqWqmppy5YtHVqSJHXVJej3A6cOTJ8CPDCqph2XPwF4uJ3/j1X1UFX9D/BJ4JyjbVqS1F2XoL8LOCvJ6UmOB7YDO1bV7AAube9fDNxRVQXcDrw0yQ+0bwA/DXx5PK1LkrrYtF5BVR1McgVNaC8AN1bV7iRXA8tVtQO4Afhokr00W/Lb2999JMkHad4sCvhkVX3iGP1bJElDpNnwnh5LS0u1vLw86TY0KzKwe2jUc3meazQ3kuyqqqVhyzwyVpJ6zqCX+sjr6WqAQS/1kdfT1QCDXuojr6erAQa91EdeT1cDDHqpj7yergYY9FIfeT1dDTDoJannDHpJ6jmDXpJ6zqCXpJ4z6DWdPLJTGhuDXtPJIzulsTHoNZ08slMaG4Ne08kjO6WxMeg1nTyyUxobg17TySM7pbEx6CWp5wx6Seo5g16Ses6gl6SeM+glqecMeknqOYNeknrOoJeknjPoJannDHppXnkq6Llh0EvzylNBzw2DXppXngp6bhj00rzyVNBzw6CX5pWngp4bBr00rzwV9Nww6CWp5zoFfZILkuxJsjfJlUOWb05yS7v8ziRbVy0/Lcl3krxrPG1LkrpaN+iTLADXARcCi8AlSRZXlV0GPFJVZwLXAtesWn4t8DdH364k6XB12aI/D9hbVStV9RhwM7BtVc024Kb2/m3A+UkCkOQNwAqwezwtS5IOR5egPxm4f2B6fztvaE1VHQQeBU5M8kzg3cD71nqAJJcnWU6yfODAga69S5I66BL0GTKvOta8D7i2qr6z1gNU1fVVtVRVS1u2bOnQkiSpq00davYDpw5MnwI8MKJmf5JNwAnAw8DLgYuTfAB4DvBEkv+tqj886s4lSZ10Cfq7gLOSnA58HdgOvGVVzQ7gUuCzwMXAHVVVwE8dKkjyXuA7hrwkbax1g76qDia5ArgdWABurKrdSa4GlqtqB3AD8NEke2m25Lcfy6YlSd2l2fCeHktLS7W8vDzpNnQsraw0Z0rcs6c5v8rOncOPyszArp9Rz1NrRteM63E0E5LsqqqlYcs8MlYbz9PjShvKoNfG8/S40oYy6LXxPD2utKEMem08T48rbSiDXhvP0+NKG8qgl6SeM+glqecMeknqOYNe0mgrK3D22bBpU3O7sjLpjnQEDHpJo3lwWy8Y9JJG8+C2XjDoJY3mwW29YNBLGs2D23rBoJc0mge39YJBL0k9Z9BLUs8Z9JLUcwa9JPWcQS9JPWfQS1LPGfSS1HMGvST1nEEvST1n0EtSzxn0ktRzBr3GywtVSFPHoNd4eaEKaeoY9BovL1Qxf/wUN/UMeo2XF6qYP36Km3oGvcbLC1XMHz/FTT2DXuPlhSrmj5/ipp5BL+no+Clu6hn0ko6On+KmnkEvST3XKeiTXJBkT5K9Sa4csnxzklva5Xcm2drO/9kku5Lc096+drztS5LWs27QJ1kArgMuBBaBS5Isriq7DHikqs4ErgWuaec/BFxUVS8BLgU+Oq7GJUnddNmiPw/YW1UrVfUYcDOwbVXNNuCm9v5twPlJUlWfr6oH2vm7gacn2TyOxiVJ3XQJ+pOB+wem97fzhtZU1UHgUeDEVTVvAj5fVd9d/QBJLk+ynGT5wIEDXXuXJHXQJegzZF4dTk2Ss2mGc9427AGq6vqqWqqqpS1btnRoSZLUVZeg3w+cOjB9CvDAqJokm4ATgIfb6VOAjwO/UlX7jrZhSdLh6RL0dwFnJTk9yfHAdmDHqpodNDtbAS4G7qiqSvIc4BPAVVX1mXE1LUnqbt2gb8fcrwBuB+4Fbq2q3UmuTvL6tuwG4MQke4F3Aoe+gnkFcCbwO0nubn+eN/Z/hSRppFStHm6frKWlpVpeXp50GzoaGdhlM+r5Zc2xr5mmXlZWmrNa7tnTnAtn506PoB2zJLuqamnYMo+MlXTseSrjiTLoJR17nsp4ogx6deeVhHSkPJXxRBn06s6P3zpSnsp4ogx6defHbx0pT2U8UQa9uvPjtzSTDHp158dvaSYZ9OrOj9/STDLoJannDHpJ6jmDXpJ6zqCXNB08IO+YMeglTQcPyDtmDHpJ08ED8o4Zg17SdPCAvGPGoJc0HTwg75gx6NVwR5gmzQPyjhmDXg13hEm9ZdCr4Y4wqbcMejXcEaZZ4BDjETHo1XBHmGaBQ4xHxKBXwx1hmgUOMR4Rg17S7HCI8YgY9JJmh0OMR8SglzQ7HGI8Igb9PPCbCponPt+fwqCfB35TQfPE5/tTGPTzwG8qaJ74fH8Kg34e+E0FzZMuz/c5G94x6OeB31TQPOnyfJ+z4R2DftZ12TLxmwqaJ12e73M2vGPQz7o52zKRxmLOhncM+mnW5Yk2Z1sm0liMa3hnRt4MOgV9kguS7EmyN8mVQ5ZvTnJLu/zOJFsHll3Vzt+T5OfG1/qM6/IE6fJEc0erdPjGNbwzrjeDY/2GUVVr/gALwD7gDOB44AvA4qqaXwM+1N7fDtzS3l9s6zcDp7d/Z2Gtxzv33HPriOzbV7W4WLWw0Nzu23d4yze6ZnGx6mlPq4LmdnHxqTULC83yQz8LC8Mfa7Bm2GNVPblmFGv6VTNNvcxizeLik2uO9DXa5bXepWYdwHKNyvFRC/6/AF4J3D4wfRVw1aqa24FXtvc3AQ8BWV07WDfq54iDfr0VNa6VPa6ark+Q9Z5oVdP14rBmemqmqZdZrOmyETWuN4MuNetYK+g3ddjoPxm4f2B6P/DyUTVVdTDJo8CJ7fx/WfW7J69+gCSXA5cDnHbaaR1aGmK9j1ldPoZtZM2LXgT33tv8tybDh1x27oRXvAIeeghOOmn01yLf857h862Z75pp6mUWa84448k1w4Z3urxGu7zWu9QcjVHvAId+gDcDfzIw/cvAH6yq2Q2cMjC9jyborwPeOjD/BuBNaz3e3GzRdxnekTT7xjXcuw7W2KLvsjN2P3DqwPQpwAOjapJsAk4AHu74u+Oxcye8+MWwsNDcDturvtbyja454wzYvRsOHmxu/W671E9dXuvHOA/SvBGsUdAE978B5wNfB+4C3lJVuwdq3gG8pKrenmQ78AtV9YtJzgY+BpwHvAD4e+Csqnp81OMtLS3V8vLyUf6zJGm+JNlVVUvDlq07Rl/NmPsVNDtSF4Abq2p3kqtpPirsoBmS+WiSvTRb8tvb392d5Fbgy8BB4B1rhbwkafzW3aLfaG7RS9LhW2uL3iNjJannDHpJ6jmDXpJ6zqCXpJ6bup2xSQ4AXz2KP3ESzSkYZsWs9Qv2vFFmredZ6xf61fMPV9WWYb8wdUF/tJIsj9rzPI1mrV+w540yaz3PWr8wPz07dCNJPWfQS1LP9THor590A4dp1voFe94os9bzrPULc9Jz78boJUlP1scteknSAINeknquN0G/3gXMp1GSryS5J8ndSabyTG5JbkzyYJIvDcx7bpJPJfn39vYHJ9njaiN6fm+Sr7fr+u4kPz/JHgclOTXJp5Pcm2R3kl9v50/tel6j52lez09P8q9JvtD2/L52/ulJ7mzX8y1Jjp90r7Bmvx9J8h8D6/hl6/6xUVckmaUfOlzAfBp/gK8AJ026j3V6fDVwDvClgXkfAK5s718JXDPpPjv0/F7gXZPubUS/zwfOae8/m+b6D4vTvJ7X6Hma13OAZ7X3jwPuBF4B3Apsb+d/CPjVSfe6Tr8fAS4+nL/Vly3684C9VbVSVY8BNwPbJtxTL1TVP9FcY2DQNuCm9v5NwBs2tKl1jOh5alXVN6rqc+39/wLupbm28tSu5zV6nlrV+E47eVz7U8Brgdva+VOzntfo97D1JeiHXcB8qp90rQL+Lsmu9gLps+KHquob0LzggedNuJ+urkjyxXZoZ2qGQQYl2Qr8OM3W20ys51U9wxSv5yQLSe4GHgQ+RTMS8K2qOtiWTFV2rO63qg6t4/e36/jaJJvX+zt9CfoMmTcL3xv9yao6B7gQeEeSV0+6oR77I+CFwMuAbwC/N9l2nirJs4C/AH6jqr496X66GNLzVK/nqnq8ql5Gc/3q84AfHVa2sV2NtrrfJD8GXAW8GPgJ4LnAu9f7O30J+o27CPkYVdUD7e2DwMdpnniz4JtJng/Q3j444X7WVVXfbF80TwB/zJSt6yTH0QTmn1XVX7azp3o9D+t52tfzIVX1LeAfaMa8n9NeGxumNDsG+r2gHTarqvou8GE6rOO+BP1dwFnt3vPjaa5Zu2PCPa0pyTOTPPvQfeB1wJfW/q2psQO4tL1/KfDXE+ylk0OB2XojU7Suk4Tmusv3VtUHBxZN7Xoe1fOUr+ctSZ7T3n8G8DM0+xY+DVzclk3Neh7R730Db/6h2Z+w7jruzZGx7de4fp/vX8D8/RNuaU1JzqDZiofmIu0fm8aek/w58BqaU6N+E3gP8Fc031Q4Dfga8OaqmpqdnyN6fg3NcELRfNvpbYfGvyctyauAfwbuAZ5oZ/82zZj3VK7nNXq+hOldzy+l2dm6QLORe2tVXd2+Fm+mGQb5PPDWdmt5otbo9w5gC82Q9d3A2wd22g7/W30JeknScH0ZupEkjWDQS1LPGfSS1HMGvST1nEEvST1n0EtSzxn0ktRz/wdq9JKteL7ypAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The number of counts is a discrete random variable.\n",
    "# Radioactive decays are randomly distributed in time, hence\n",
    "# the probability distribution for the number of counts is a Poisson with mean 18.0\n",
    "counts = stats.poisson(xbar)\n",
    "\n",
    "# plot the pmf\n",
    "x = np.arange(35)\n",
    "plt.plot(x,counts.pmf(x), 'ro', ms=4)\n",
    "plt.vlines(x, 0, counts.pmf(x), colors='r', lw=3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) What is the probability of observing fewer than 10 counts within a ten-minute interval?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.015381097260589295"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The probability of observing 9 or fewer counts is given by the CDF:\n",
    "counts.cdf(9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "The G-M tube has a circular window with diameter 5cm.\n",
    "\n",
    "The tube is positioned 10cm from the centre of the charcoal sample.\n",
    "\n",
    "The tube’s detection efficiency for the low-energy beta particles emitted by C-14 is 2%.\n",
    "\n",
    "---\n",
    "\n",
    "### Question 3 \n",
    "\n",
    "Calculate a 90% confidence interval for the mean number of beta decays occurring in the sample per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "window area:  0.001963495408493621\n",
      "sphere surface area:  0.12566370614359174\n",
      "fraction of decays detected:  0.0003125\n",
      "mean decays/s: 96.0\n",
      "stderr decays/s: 11.845924443703241\n"
     ]
    }
   ],
   "source": [
    "# area of the window:\n",
    "a = np.pi * (0.025)**2\n",
    "print('window area: ', a)\n",
    "\n",
    "# surface area of a 10cm sphere around the sample:\n",
    "A = 4 * np.pi * (0.1)**2\n",
    "print('sphere surface area: ', A)\n",
    "\n",
    "# efficiency of the detector:\n",
    "eff = 0.02\n",
    "\n",
    "# fraction of decays that are actually detected:\n",
    "f = (a/A) * eff\n",
    "print('fraction of decays detected: ', f)\n",
    "\n",
    "# converting the mean and standard error to total decays per second:\n",
    "mean_decays = xbar / ( 10 * 60 * f )\n",
    "print('mean decays/s:', mean_decays)\n",
    "\n",
    "stderr_decays = stderr / ( 10 * 60 * f )\n",
    "print('stderr decays/s:', stderr_decays)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% CI for decays per second: [ 72.13 119.87]\n"
     ]
    }
   ],
   "source": [
    "# Because n<30, we need to use Student's t-distribution to find the confidence interval.\n",
    "T = stats.t(df=n-1)\n",
    "t90 = T.interval(0.9)[1]\n",
    "\n",
    "ci = np.array([mean_decays - t90 * stderr_decays , mean_decays + t90 * stderr_decays])\n",
    "print('90% CI for decays per second:',ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "The charcoal is estimated to be 80% carbon by weight.\n",
    "\n",
    "A carbon atom has a mass of $2 \\times 10^{-26}$ kg.\n",
    "\n",
    "Approximately 1 in every $10^{12}$ carbon atoms in the atmosphere is C-14.\n",
    "\n",
    "The [half-life](https://en.wikipedia.org/wiki/Half-life) of a radioisotope is the period of time after which its probability of having decayed is 0.5.\n",
    "\n",
    "The half-life, $t_{1/2}$, is given by the formula\n",
    "$$ t_{1/2} = \\tau\\ln (2), $$\n",
    "where $\\tau$ is the *mean lifetime* of the radioisotope.\n",
    "\n",
    "---\n",
    "\n",
    "### Question 4\n",
    "\n",
    "Using your answer to Q3, calculate a 90% confidence interval for the half-life of Carbon-14."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sample contains 2.e+13 atoms of C-14\n",
      "90% CI for decays per C-14 atom per second: [3.61e-12 5.99e-12]\n"
     ]
    }
   ],
   "source": [
    "# First, we should convert our confidence interval into decays per C-14 atom per second:\n",
    "\n",
    "carbon_mass = 0.5 * 0.8\n",
    "n_C = carbon_mass / 2e-26\n",
    "n_C14 = n_C / 1e12\n",
    "print('sample contains', np.format_float_scientific(n_C14), 'atoms of C-14')\n",
    "\n",
    "\n",
    "ci_C14 = ci / n_C14\n",
    "print('90% CI for decays per C-14 atom per second:',ci_C14)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% CI for mean lifetime (in seconds): [1.67e+11 2.77e+11]\n",
      "90% CI for C-14 half-life (in seconds): [1.16e+11 1.92e+11]\n"
     ]
    }
   ],
   "source": [
    "# Now we can use the formula given to calculate the CI for half life.\n",
    "\n",
    "# The mean lifetime is the inverse of the decay constant (expected decays per second, aka lambda).\n",
    "# We can plug in the confidence interval that we just calculated.\n",
    "# The numpy flip() function reverses the order of the CI array.\n",
    "tau = np.flip(1 / ci_C14)\n",
    "print('90% CI for mean lifetime (in seconds):', tau)\n",
    "\n",
    "t_half_ci = tau * np.log(2)\n",
    "print('90% CI for C-14 half-life (in seconds):', t_half_ci )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% CI for C-14 half-life (in years): [3664.72 6090.27]\n"
     ]
    }
   ],
   "source": [
    "# or in years:\n",
    "print('90% CI for C-14 half-life (in years):', t_half_ci/(60*60*24*365.25) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "The sample is replaced by a new charcoal sample, of the same mass but unknown origin. \n",
    "\n",
    "The G-M counter reads the following counts over six consecutive ten-minute periods: 8, 16, 10, 20, 17, 14.\n",
    "\n",
    "---\n",
    "\n",
    "### Question 5\n",
    "\n",
    "Is this sample older than the first one? Carry out a hypothesis test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Difference in mean counts: 3.833333333333334\n",
      "p-value: 0.10636637637680363\n"
     ]
    }
   ],
   "source": [
    "data2 = np.array([8, 16, 10, 20, 17, 14])\n",
    "\n",
    "# If the second sample is substantially older, it will have lower radioactivity as some of its C-14 will\n",
    "# have already decayed.\n",
    "\n",
    "# We can use an independent 2-sample t-test to investigate the difference in the mean number of counts,\n",
    "# d = mu_sample1 - mu_sample2\n",
    "# H0 : d = 0\n",
    "# H1 : d > 1\n",
    "\n",
    "# Let's choose a significance level alpha = 0.05\n",
    "\n",
    "# The difference in means is\n",
    "d = data.mean() - data2.mean()\n",
    "print('Difference in mean counts:',d)\n",
    "\n",
    "# In the course notebook, we did the t-test the long-winded way.\n",
    "#\n",
    "# Actually, scipy has a convenient t-test function, but we need to remember that it returns a \n",
    "# TWO-TAILED p-value. Always check the documentation!\n",
    "# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html\n",
    "\n",
    "# We feed in the two data samples and collect the two-tailed p-value:\n",
    "p_2tailed = stats.ttest_ind(data,data2).pvalue\n",
    "\n",
    "# then halve the result to give the one-tailed p-value:\n",
    "p = p_2tailed / 2\n",
    "print('p-value:',p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# So there is no significant difference in the sample means, at the 5% significance level.\n",
    "# Hence we cannot reject H0. There is no evidence that the second sample is older than the first."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
